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1.  Introduction 


Sun  (2011)  introduced  a  new  semi-implieit,  time  integration  scheme  and  simulated 
solutions  to  the  eylindrical  dam-break  problem  for  the  nonlinear  shallow  water 
equations.  The  seheme  exhibited  several  important  benefits  relating  to  the  overall 
stability  of  the  integration  and  preservation  of  strong  gradients  in  the  calculated 
flow  quantities.  The  integration  reproduced  the  solutions  of  Godunov-type 
Riemann  solvers  eomparing  favorably  to  more  sophisticated  and  more 
computationally  intensive  schemes.  The  new  scheme  is  a  2-time-step  scheme 
similar  in  structure  to  the  forward-backward  finite-difference  seheme  (Haltiner  and 
Williams  1980).  Thus,  there  is  no  eomputational  mode,  and  ineorporating  a  split- 
time-step  seheme  in  order  to  less  frequently  update  slowly  varying  quantities  will 
not  require  any  time  filtering  (Hsu  and  Sun  2001).  As  part  of  the  linear  stability 
analysis.  Sun  (2011)  showed  that  the  scheme  has  a  diffusion  term  similar  to  the 
Lax-Wendroff  seheme;  however,  based  on  eigenvalue  analysis,  the  amplification 
factor  is  unity  for  Courant  numbers  less  than  1 . 

Here  we  apply  the  seheme  to  the  3-dimensional  (3-D)  Navier-Stokes  equations  as 
a  module  in  the  Atmospherie  Boundary  Layer  Environment  (ABLE)  model  (Wang 
et  al.  2012),  a  new  eomputational  fluid  dynamics  (CED)  model  eode  in 
development  at  the  US  Army  Research  Eaboratory.  The  model  is  being  developed 
for  eventual  applieation  to  atmospheric,  planetary  boundary  layer  (PBE)  flows  in 
complex  and  urban  terrain.  Aeeurate  simulation  within  the  PBE  necessarily 
includes  producing  fine-scale  struetures  (e.g.,  thermal  plumes  in  eonvective 
conditions  [Moeng  and  Sullivan  1994],  shear-layers  [Chimonas  1999],  andterrain- 
indueed  vortex  shedding/generation  in  stable  eonditions  [Grubisic  et  al.  2008]),  and 
exeessive  numerieal  diffusion  may  eompromise  their  formation  or  severely  degrade 
their  evolution  compared  to  observations. 

The  new  semi-implieit  seheme  was  implemented  as  an  option  for  the  ABLE 
model’s  dynamieal  eore  (version  0.90).  In  order  to  characterize  the  scheme’s 
performanee,  we  simulated  a  3-D  lid-driven  cavity  flow  (reviewed  in  Shankar  and 
Deshpande  [2000])  and  eompared  the  results  simulations  of  Ku  et  al.  (1987)  and 
the  laboratory  measurements  of  Prasad  and  Koseff  (1989).  Eid-driven  cavity  flow 
is  commonly  used  to  evaluate  numerical  schemes  because  of  the  unambiguous 
boundary  eonditions,  the  detailed  laboratory  observations  of  Prasad  and  Koseff 
(1989),  and  features  such  as  the  intermittent  and  semi-permanent  vortex  struetures 
shown  in  Pig.  1  (Shankar  and  Deshpande  2000). 
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Fig.  1  Diagram  of  lid-driven  cavity  flow  and  associated  flow  structnres.  The  npper 
bonndary  of  the  hox  is  moving  with  a  velocity  of  I/q  in  the  positive  x-direction.  All  honndaries 
are  no-slip.  Adapted  from  Shankar  and  Deshpande  (2000). 


2.  Methodology 


The  general  formulation  for  the  semi-implicit  scheme  for  a  set  of  n  prognostic 
variables,  0j,  impacted  by  forcing  functions,  7^,  is 


dt 


i  =  1,  ...,n;  j  =  1,2,3. 


(1) 


Discretizing  the  above  set  of  n  equations,  the  semi-implicit  scheme  calculates  the 
forcing  at  the  midpoint  between  the  2  time  steps,  t  =  to  +  At/2; 


4>i{xjAx.+to)-4>i{xj,to) 

At 

1,  ...,n;  i  =  1,2,3 


(2) 


Approximating  the  prognostic  variables  at  that  midpoint  is  done  using  a  truncated 
Taylor  expansion  with  respect  to  time.  The  time  derivatives  in  the  first-order  terms 
are  then  replaced  by  the  forcing  functions,  ^Fj,  evaluated  at  t  =  tg: 


/  At  \  ,  .  Atd0j(xy,t) 

0t  ~  +y - ^ - 


t  =  to 


(3) 


For  the  3-D  finite-volume  formulation,  the  prognostic  variables  are  density,  p, 
momentum,  l/j  =  pui,  and  density-coupled  potential  temperature,  Q  =  p6.  For  the 
simulations  reported  below,  the  model  is  configured  as  a  direct  numerical 
simulation  (DNS).  Because  the  molecular  diffusion  terms  are  retained,  runs  must 
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be  high  enough  resolution  to  resolve  motions  down  into  the  dissipation  subrange, 
where  the  magnitudes  of  these  diffusion  terms  beeome  of  similar  order  to  the  other 
foreing.  Beeause  of  the  importanee  of  aeeurately  ealeulating  the  pressure  gradient 
foree  in  the  momentum  equations,  the  thermodynamie  variables,  p  and  0,  are 
updated  first.  Following  the  summation  eonvention  over  repeated  indiees,  the 
governing  equations  are 


where  the  pressure  is  ealeulated  diagnostieally  from  the  Ideal  Gas  Law 
P^vi^v=RQ/p^  p  ^  slows  the  propagation  speed  of  aeoustie  waves 

(c  =  ^JyWtJS'),  allowing  for  longer  time  steps  (Sun  et  al.  2012).  R  is  the  Universal 
Gas  Constant  divided  by  the  mean  molar  mass  of  the  fluid,  and  Cp  and  Cy  are  the 
appropriate  speeifie  heats  at  eonstant  pressure  and  volume,  respeetively.  The  other 
variables  are  p,  the  moleeular  viseosity;  k,  the  moleeular  eonduetivity;  and  pi,  the 
veetor  gravitational  aeeeleration. 

Using  the  prognostie  density  equation  requires  temporally  resolving  the  aeoustie 
waves;  thus,  simulations  employ  exeeedingly  small  time  steps  for  integration 
the  physieal  time  seale).  The  simulations  deseribed  in  this  study 
employed  S  =  16,  thus  stable  integration  is  aehieved  using  a  time  step  4  times 
larger  than  normally  required  without  signifieantly  impaeting  flow  struetures  larger 
than  a  few  times  the  grid  spaeing  (Sun  et  al.  2012).  A  time-splitting  seheme,  where 
the  more  slowly  varying  foreing  meehanisms  are  updated  less  often  (e.g.,  Hsu  and 
Sun  2001),  is  another  alternative  to  deereasing  the  required  eomputational 
resourees. 

The  new  model  uses  the  quadratie  upstream  interpolation  for  eonveetion  kineties 
(QUICK)  finite-volume  eonveetion-diffusion  seheme  (Versteeg  and  Malalasekera 
2007)  modified  to  allow  non-uniform  grid  spaeing.  The  volume  elements  for  the 
eomponents  of  the  momentum  field  are  staggered  with  respeet  to  the  mass  field 
(Arakawa  C-grid).  The  volume  elements  are  eonfigured  so  that  eomplete  mass-field 
volume  elements  are  in  eontaet  with  the  wall  (Fig.  2).  The  outermost  staggered 
veloeity  volume  elements  are  between  the  2  outermost  mass-field  eells  in  the 
direetion  of  the  veloeity.  In  this  eonfiguration,  no  pressure  boundary  eondition  is 
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required  for  the  momentum  equation.  Diriehlet  boundary  eonditions  are  imposed 
on  all  prognostie  fields. 


•  Scalar-cell  center 

•  U-cell  center 

•  V-ceU  center 


Fig.  2  Arakawa  C-grid  configuration  in  the  xy-plane  for  the  scalar  mass  field  and  the  U- 
and  V-components  of  the  momentum  field.  The  vertical  momentum  is  staggered  in  the  vertical 
direction  from  the  mass  field  points. 

Using  the  above  semi-implieit  seheme  yields  the  following  governing  equations  in 
Cartesian  eoordinates  with  fluxes  eomputed  on  the  surfaee  of  finite  volume 
element: 
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p(x,At  +  to)  -p(x,to)  =  -  ASj-^1/^-^) 


(V) 


0(x,At+  to)  -  0(x,to)  = 


At 

AV 


i:, 


3(-) 


dxj  p(  ) 


(8) 


At 


0j(x,At  +  to)  -  0i(x,to)  =  — 


I', 


as; 


(_)  f  a  <£/;>(-)' 


K-) 


dxj  p(  ) 


+ 


(+)/  a  <Ui)W 


As;-M/t 


dxj  pW 


)(+) 


+  Mpgt  -  ^(pWAS«  -  P«ASJ-^) 


(9) 


where  the  indiees,  zj=l,2,3,  indicate  the  vector  components  and  cell  faces  in  each 
direction.  The  superscripts  (+)  and  (-)  indicate  values  on  the  cell  faces  in  the 
positive  or  negative  direction  from  the  cell  center;  AS[^^  and  AS[  ^  are  the  surface 
elements  of  the  cell  faces  in  the  /-direction.  With  the  staggered  grid,  the  momentum 
flux  contributions  in  the  continuity  equation  and  the  pressure  gradient  term  in  the 
momentum  equation  are  directly  calculated  rather  than  using  interpolated  values. 
Two  types  of  interpolation  are  used  for  the  convection-diffusion  terms.  An  over  bar 
indicates  simple  linear  averaging  between  grid  cell  centers,  if  necessary;  because 
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of  the  cell  staggering  employed,  not  all  faces  will  require  interpolation  to  find 
values  on  the  faces.  The  ( )  operator  indicates  quadratic,  upwind-weighted 
interpolation  between  surrounding  cell  centers  in  the  direction  of  the  flow  through 
the  face,  as  discussed  in  Versteeg  and  Malalasekera  (2007). 


The  integration  process  is  described  using  the  above  equations; 


1) 


Using  the  discretized  momentum  equation  (Eq.  9),  calculate  — 
estimate  the  half-step  momentum,  l/j(x,  tg  +  y)- 


t=to 


in  order  to 


2)  With  the  half-step  momentum,  advance  the  density  a  full  time  step  (Eq.  7),  and 
calculate  the  half-step  density  as  a  linear  average  between  the  old  and  updated 
densities. 


3)  The  density-coupled  potential  temperature  is  then  advanced  by  a  half  time  step 
using  the  first-order  Taylor  expansion;  the  full-time-step  advancement  follows 
employing  the  updated  mid-time-step  density  and  density-coupled  potential 
temperature. 

4)  With  the  updated  scalar  quantities,  the  mid-time-step  pressure  is  calculated,  and 
the  momentum  equation  is  advanced  to  the  next  full  time  step. 

5)  Einally,  as  recommended  in  Sun  (2011),  a  light  application  of  fourth-order 
Shuman  smoothing  (e.g.,  Haltiner  and  Williams  [1980]  or  Shapiro  [1975]) 
eliminates  short-wave  instability. 


We  found  stable  time  integration  for  all  simulations  using  20%  of  the  maximum 
smoothing  coefficient  and  smoothing  every  16  time  steps. 


3.  Results 


The  ABEE  model,  using  the  above  governing  equations,  simulated  3-D  lid-driven 
cavity  flow;  results  are  compared  with  the  laboratory  measurements  of  Prasad  and 
Koseff  (1989).  We  specifically  looked  at  the  cases  where  the  stream- wise,  span- 
wise  and  vertical  lengths  were  all  equal  (SOR=l)  and  Reynolds  numbers  were 
7?e=1000,  7?e=3200,  and  7?e=10,000.  The  results  are  normalized  by  the  flow  scales; 
the  length  scale,  L,  is  the  cube  width;  the  velocity  scale,  Uq,  is  the  speed  of  the 
moving,  upper  lid;  and  the  time  scale  is  t  =  L/Uq. 


3.1  Discretization  Error 

The  theoretical  accuracy  of  the  numerical  schemes  has  been  derived  in  Sun  (2011) 
and  Versteeg  and  Malalasekera  (2007)  as  second-order  accurate  in  time  and,  when 
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using  the  QUICK  scheme,  third-order  accurate  in  space.  As  implemented  in  the 
model  source  code,  the  spatial  accuracy  decreases  due  to  the  boundary  condition 
treatment  in  both  the  convection-diffusion  scheme  and  the  smoothing  routine. 

The  temporal  accuracy  was  measured  by  simulating  cavity  flow  for  i?e=  10,000 
using  a  mesh  with  304x304x304  grid  points.  The  simulation  used  5=16  for  a  4 
times  reduction  in  the  speed  of  sound;  the  largest  time  step  satisfied  the  Courant- 
Freidrichs-Levy  (CFL)  criterion  for  the  slowed  acoustic  waves  with  a  value  of 
5.8x1 0'^r.  Three  additional  simulations  were  performed  using  time  steps  of  1/2, 1/4, 
and  1/8  of  the  largest  time  step.  The  T2-norm  of  the  difference  between  each 
simulation  and  the  simulation  with  the  smallest  time  step  taken  at  t=5.ST  (10^ 
iterations  of  the  largest  timestep)  is  shown  in  Fig.  3  along  with  the  line  representing 
second-order  accuracy  in  time. 


Fig.  3  Estimate  of  temporal  discretization  error  based  on  lid-driven  cavity  flow  simulations 
with  f2e=10,000.  The  largest  time  step  satisfied  the  CFL  criterion  for  the  reduced  speed 
acoustic  waves.  Additional  simulations  were  performed  using  an  identical  grid,  but  with  time 
steps  reduced  by  factors  of  1/2, 1/4,  and  1/8,  which  was  used  as  the  reference  state. 

For  the  spatial  error  estimates,  simulations  were  performed  with  1000  to  reduce 
the  computational  requirements.  The  QUICK  scheme  is  third-order  accurate  in 
space.  However,  the  boundary  condition  treatment  for  both  the  convection- 
diffusion  scheme  and  the  explicit  smoothing  degrades  the  accuracy.  Four 
simulations  with  differing  grid  spacings  were  performed;  75, 150,  300,  and  450  grid 
points  were  used  in  each  Cartesian  direction.  The  time  steps  for  all  of  the 
simulations  were  set  at  the  time  step  satisfying  the  CFL  criterion  for  the  450-grid 
point  mesh,  and  the  simulations  were  allowed  to  integrate  t=40T  (~10^  time  steps). 
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The  coarser  resolution  simulations  were  compared  with  the  finest  resolution  using 
third-order  splines  to  interpolate  between  grid  points.  Mean  absolute  error  of  the 
difference  between  the  interpolated  values  and  the  450-grid-point  simulation  were 
then  calculated  over  the  area  at  the  same  physical  time,  t=40r,  and  plotted  in  Fig.  4. 


Fig.  4  Spatial  discretization  error  estimates  using  stream-wise  and  vertical  velocities  in  the 
plane  of  symmetry  (y=0)  for  f?e=1000.  The  error  is  the  Ti-norm  of  the  variable  difference  from 
a  reference  simulation  that  used  450  grid  points  in  each  direction.  For  a)  the  error  is  calculated 
over  the  entire  symmetry  plane  and  for  h)  20yo  of  each  side  of  the  symmetry  plane  is  excluded. 

The  error  estimates  are  calculated  for  the  stream-wise  and  vertical  velocities  for  2 
regions  of  the  central  symmetry  plane  (y=0).  The  first  region  is  the  complete  xz- 
slice  through  the  center  of  the  computational  domain.  The  second  region  excludes 
20%  of  the  grid  points  along  each  boundary.  Comparing  Fig.  4a  and  4b,  it  is  readily 
apparent  that  most  of  the  error  results  from  the  boundary  condition  treatment;  the 
larger  magnitude  of  the  upper  boundary  condition  of  the  w-velocity  contributes 
more  to  the  mean  absolute  error.  The  QUICK  scheme  is  a  third-order  scheme,  but 
the  second-order  treatments  for  both  convection-diffusion  and  second-order 
smoothing  at  the  wall  boundaries  degrades  the  discretization  error  to  match  second 
order  in  space. 
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3.2  Comparison  with  Simulation  and  Observations 

The  simulations  were  conducted  as  DNSs  requiring  high  resolution  in  order  to 
resolve  motions  into  the  dissipation  sub-range:  at  least  75,  200,  and  300  grid  points 
in  each  direction  for  7?e=  1000, 7?e=3200,  and  7?e=  10,000,  respectively.  Gravity  was 
disabled  for  these  simulations  and  the  potential  temperature  was  fixed. 

For  i?e=1000,  the  new  model  simulations  are  compared  with  the  simulation  data  of 
Ku,  Hirsh,  and  Taylor  (1987),  which  employed  a  spectral  DNS  to  solve  for  the 
3-D  flow  field.  Their  3-D  simulations  showed  a  significant  change  in  the  profile 
when  compared  with  the  horizontal  and  vertical  centerline  velocity  profiles 
generated  using  2-dimensional  (2-D)  simulations.  Fig.  5  shows  the  stream-wise 
velocity  along  the  central  vertical  axis  of  the  symmetry  plane  and  the  vertical 
velocity  along  the  central  horizontal  axis.  All  the  simulations  were  in  close 
agreement  with  the  simulations  of  Ku,  Hirsh,  and  Taylor  (1987). 
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(b)  X* 

Fig.  5  a)  Stream-wise  velocity  along  the  central  vertical  axis  and  b)  vertical  velocity  along 
the  central  horizontal  axis  for  i?e=1000.  The  model  simnlations  nsed  75,  150, 300,  and  450  grid 
points  in  each  direction  and  are  compared  with  the  simnlations  of  Kn,  Hirsh,  and  Taylor 
(1987). 

The  simulations  at  higher  Reynolds  number  (Re  =  3200  and  Re  =  10,000)  are 
verified  against  the  laboratory  observations  of  Prasad  and  Koseff  (1989).  For  better 
comparison  with  the  statistics  of  the  Doppler  velocimetry  observations,  the  points 
along  the  central  vertical  axis  in  the  symmetry  plane  (x  =  y  =  0)  and  the  central 
horizontal  axis  (y  =  z  =  0)  in  the  symmetry  plane  were  written  for  every  time  step. 
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which  was  ~10  (i?e=3200)  and  ~5  x  10  (Re=10,000).  For  both  cases,  the 

model  was  integrated  for  a  physieal  time  of  approximately  120t. 

A  representative  sample  of  the  evolution  of  the  flow  for  i?e=3200  is  shown  in  the 
density  and  streamline  plots  of  Fig.  6.  The  simulation  was  allowed  to  aehieve  quasi- 
steady-state,  which  exhibited  a  periodie  shedding  of  the  downstream  seeondary 
eddy  (DSE)  and  the  upstream  seeondary  eddy  (USE)  (see  Eig.  1).  The  flow  shows 
the  prototypieal  state  with  a  large  primary  eddy  (PE)  and  well-defined  seeondary 
eddies — the  DSE,  the  USE,  and  a  small  deviation  for  the  upper  upstream  eddy 
(UUE)  (see  Eig.  1).  The  DSE  is  then  shed  into  the  primary  flow,  eausing  the 
downstream  flow  to  wrap  around  the  PE  further  from  the  lower  wall.  As  the  DSE 
is  absorbed  into  the  PE  and  USE  is  shed  into  the  primary  cireulation,  the  upward 
flow  near  the  bottom  boundary  weakens,  leading  to  a  rapid  drop  in  the  level  that 
the  downstream  flow  beeomes  horizontal.  This  drop  is  assoeiated  with  an 
aeeeleration  of  the  downstream  flow,  which  reestablishes  the  DSE  and  then  the 
USE.  Eor  the  higher-i?e  flow,  there  is  a  similar  DSE  and  USE  shedding  proeess;  the 
spike  in  the  vertieal  veloeity  variance  is  lower  magnitude  and  narrower  due  to  the 
deereased  viscous  diffusion. 

The  flow  took  approximately  35t  to  aehieve  quasi-steady-state,  which  forms  the 
lower  boundary  of  the  window  used  to  ealeulate  the  mean  and  seeond-order 
statistics.  Comparing  the  simulation  results  to  the  observations,  we  see  exeellent 
agreement  with  the  mean  features  (Eig.  7).  The  sharp  gradients  in  the  mean  field, 
espeeially  near  the  walls,  are  well  represented.  The  downstream  flow  (near  the  right 
wall)  is  highly  loealized  in  both  simulations,  with  the  expeeted  wider  profile  for  the 
lower-Re  flow. 
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Fig.  6  Streamlines  and  UW-velocity  magnitudes  (^e=3200)  showing  a)  at  t  =  90t,  a  well- 
defined  DSE  and  localized  downstream  flow  along  the  right,  boundary;  h)  at  t  =  93t,  the 
shedding  of  the  DSE  raising  the  level  that  the  downstream  maximum  wraps  around  the  PE; 
c)  at  t  =  97t,  drop  and  acceleration  of  the  downstream  flow;  and  d)  at  t  =  104t, 
reestablishment  of  the  DSE  and  thinning  of  the  downstream  maximum. 
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Fig.  7  Non-dimensional  mean  horizontal  and  vertical  velocity  Helds  along  the  x-  and  z-axes 
in  the  symmetry  plane  (y=0)  for  a)  f2e=3200  and  h)  f2e=l  0,000 

The  features  of  the  velocity  variances  (Fig.  8,  scaled  for  clarity)  also  capture  the 
observed  behavior  described  by  Prasad  and  Koseff  (1989).  The  high-peaks  in  the 
vertical  velocity  variance  near  the  right  lateral  boundary  are  associated  with  the 
thickening  and  thinning  of  the  high-velocity  shear  layer  associated  with  the 
shedding  of  the  DSE  vortex  in  the  comer  described  above.  Using  a  lower  order 
convection-diffusion  scheme  (such  as  upwind-differencing)  with  excessive 
numerical  diffusion  did  generate  the  secondary  vortices;  however,  the  shedding 
process  was  eliminated.  The  flow,  instead,  generated  a  stable  PE  with  steady  DSE 
and  USE.  The  shedding  of  the  DSE  (and  USE)  has  impacts  as  discussed  in  Prasad 
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and  Koseff  (1989):  1)  sharp  increase  in  the  veloeity  varianee  along  the  lower 
boundary  due  to  said  eddies  erossing  the  eenterline  and  2)  the  thiekening  and 
thinning  of  the  high-veloeity  downward  boundary  layer  along  the  vertieal  wall  on 
the  downstream  side,  appearing  as  a  spike  in  the  vertieal  veloeity  varianee.  The 
greater  turbulent  diffusion  assoeiated  with  the  higher-Re  flow  leads  to  a  less 
loealized  inerease  in  the  horizontal  veloeity  varianee  along  the  bottom  boundary. 


0.4  OJ,  0-0  0J2  0.4 


(b)  ^ .  8£„/f/o 

Fig.  8  Scaled  horizontal  and  vertical  velocity  variances  along  the  x-  and  z-axes  in  the 
symmetry  plane  (y=0)  for  a)  t?e=3200  and  h)  0,000 

The  effeet  of  the  shedding  of  the  DSE  and  USE  also  has  signifieant  effeet  on  the 
veloeity  covariances  (Eig.  9,  sealed  for  elarity)  produeing  a  large  positive  value 
along  the  right  lateral  boundary.  An  inerease  in  downward  veloeity  (-w’)  is 
assoeiated  with  the  thiekening  (in  the  negative  x-direetion)  for  the  downstream  flow 
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(-u’).  Similarly,  the  slowing  of  the  flow  (+w’)  is  associated  with  the  rightward 
motion,  thinning  the  downstream  flow  (+u’). 


(a) 


(b) 


Fig.  9  Scaled  horizontal  and  vertical  velocity  covariances  along  the  central  axes  in  the 
symmetry  plane  (y=0)  for  a)  f?e=3200  and  h)  f?e=10,000 

The  large  positive  peak  in  the  uw-covariance  near  the  bottom  boundary  for  the  low¬ 
ing  simulation  was  described  in  Prasad  and  Koseff  (1989)  as  being  caused  by  the 
meandering  of  the  Taylor-Goertler-like  (TGL)  vortices  on  the  Re=3200  flow 
(Fig.  10).  As  the  TGL  vortices  approach  the  symmetry  plane,  the  vertical  velocity 
increases,  and  the  horizontal  velocity  maximum  moves  upwards,  yielding  a 
decrease  in  the  negative,  x-component  of  the  velocity  (i.e.,  +u’)  relative  to  the  mean 
near  the  surface  boundary,  and  an  increase  in  the  negative  U-velocity  above 
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(i.e.,  -u’).  Coupled  with  the  upward  motion  (+w’),  the  movement  of  the  TGL 
vortiees  close  to  the  symmetry  plane  contributes  to  a  sharp  positive  increase  in  the 
uw-covariance  below  the  mean  U-velocity  maximum.  Above  the  mean  U-velocity 
maximum,  the  covariance  will  exhibit  a  negative  peak.  The  observations  from 
Prasad  and  Koseff  (1989)  show  a  meandering  period  of  about  70t,  whereas  the 
simulations  showed  a  much  smaller  meandering  period  of  approximately  14t. 

(a)  (b) 


r*  r* 


Fig.  10  Streamline  and  VW-velocity  magnitudes  for  the  same  times  shown  in  Fig.  6 
a)  t  =  90t,  h)  t  =  93t,  c)  t  =  97t,  and  d)  t  =  104t.  This  progression  of  figures  shows  the 
meandering  of  the  TGL  vortices  near  the  symmetry  plane  (Y*=0). 

For  the  i?e=10,000  flow,  the  TGL  vortices  have  a  less  regular  configuration  and  are 
transient  (Fig.  1 1).  In  addition,  the  mean  U-velocity  maximum  is  closer  to  the  lower 
boundary.  There  is  a  small  positive  incursion  in  the  uw-covariance  near  the  lower 
boundary;  however,  during  DSE  shedding  events  most  of  the  lower  part  of  the 
domain  experiences  negative  x-acceleration  (increase  in  negative  U-velocity),  with 
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an  increase  in  the  upward  component  or  a  positive  x-aeceleration  (deerease  in 
negative  U-velocity)  during  the  downward  shift  of  the  veloeity  maximum. 


Fig.  11  Streamlines  and  VW-velocity  magnitude  for  Re=10,000  at  time,  t  —  IIOt.  For  the 
higher-Re  flow  the  vortices  generated  near  the  lower  boundary  are  irregular  and  transient. 

4.  Conclusion 


We  extended  the  time  integration  scheme  of  Sun  (2011)  to  3-D  Navier-Stokes  and 
implemented  the  seheme  in  the  ABLE  model.  The  scheme  was  straightforward  to 
implement  and  parallelize  using  Message  Passing  Interface  (MPI)  (optimizing 
ABLE  for  heterogeneous  platforms  (e.g.,  general  purpose  computing  on  a  graphics 
processing  unit  [GPGPU]  is  in  progress).  The  model  showed  favorable  stability  and 
fidelity  eharaeteristies  when  applied  to  lid-driven  cavity  flow  at  i?e=1000, 
i?e=3200,  and  10,000.  The  simulations  accurately  reproduced  the  mean  velocity 
fields  and  the  second-order  statisties  along  the  symmetry  plane  (y=0)  when 
eompared  with  other  3-D  simulations  of  Ku,  Hirsh,  and  Taylor  (1987)  and  the 
laboratory  observations  of  Prasad  and  Koseff  (1989).  The  initially  imposed,  random 
perturbations  in  the  density  field  (initially,  ~10“^po)  maintained  beeause  of  the 
lower  numerical  diffusion,  despite  the  smooth  nature  of  the  boundary  eonditions, 
the  use  of  Shuman  smoothing  to  control  short-wave  instability,  and  the  small  time 
steps  required.  Maintaining  small  random  perturbations  in  the  flow  preserved  the 
transient  nature  of  the  flow  struetures,  sueh  as  the  formation  and  shedding  of  the 
DSE,  USE,  and  TGL  vortices.  Associated  with  the  shedding  of  the  corner  eddies  is 
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the  alternating  thiekening  and  thinning  of  the  lateral  boundary  layer,  leading  to  the 
spikes  in  the  vertieal  veloeity  varianee,  (T^^,  above  the  DSE  and  above  the  lower 
boundary.  The  production  and  meandering  of  TGL  vortices  in  the  Re=3200  flow 
leads  to  a  large  positive  value  for  the  uw-covariance,  as  observed  in  the  laboratory. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


2-D 

2-dimensional 

3-D 

three-dimensional 

ABLE 

Atmospheric  Boundary  Layer  Environment 

CFD 

computational  fluid  dynamics 

CEL 

Courant-Ereidrichs-Eevy 

DNS 

direct  numerical  simulation 

DSE 

downstream  secondary  eddy 

DSRC 

Defense  Shared  Resource  Centers 

GPGPU 

general  purpose  computing  on  a  graphics  processing  unit 

HPC 

high-performance  computing 

HPCMO 

High-Performance  Computing  Modernization  Office 

MPI 

Message  Passing  Interface 

PBE 

planetary  boundary  layer 

PE 

primary  eddy 

QUICK 

quadratic  upstream  interpolation  for  convective  kinetics 

TGL 

Taylor-Goertler  like 

USE 

upstream  secondary  eddy 

UUE 

upstream  upper  eddy 
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